II. ALPHA DIVERSITY ANALYSES

Loading libraries

library(hilldiv)
library(ggplot2)
library(ggpubr)
library(reshape2)
library(tidyverse)
library(lme4)
library(car)
library(arm)
library(predictmeans)
library(AER)
library(visreg)
library(FSA)
library(rstatix)
library(sjPlot)
library(effects)
library(cowplot)

Calculating Hill numbers with hilldiv R package

otu_table <- read.csv("feature_table.csv", header = TRUE, row.names = 1)
q0 <- hill_div(otu_table, qvalue = 0)
q1 <- hill_div(otu_table, qvalue = 1)
q2 <- hill_div(otu_table, qvalue = 2)
q012 <- cbind(q0, q1, q2)
write.table(q012, file="../data/q012_hilldiv.txt", sep = "\t")

Preparing the data

Load files

richness_q012 <- read.csv("../data/Hill_numbers_q012.csv", header = TRUE) %>% 
  dplyr::select(SampleID, q0, q1, q2)
metadata <- read.csv("../data/metadata.csv", header = TRUE, check.names = F)
Micro_div <- richness_q012 %>% 
  inner_join(metadata, by = c("SampleID"="SampleID"))

Declare factors

SPECIES <- as.factor(Micro_div$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(Micro_div$Season) # Dry and Rainy
SEX <- as.factor(Micro_div$Sex) # Male and Female

Standardize continous independent variables

SVL <- rescale(Micro_div$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(Micro_div$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(Micro_div$SeqDepth, binary.inputs = "center")

Evaluate whether sex variable influence on gut microbiota diversity (q=1)

# Paired test (Wilcoxon test, q1) 

# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                            Sex == "Female", q1, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                             Sex == "Male", q1, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  aeneus_fem_GutMic and aeneus_male_GutMic
## W = 49, p-value = 0.4779
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                          Sex == "Female", q1, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                           Sex == "Male", q1, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                          Sex == "Female", q1, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                           Sex == "Male", q1, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  gram_fem_GutMic and gram_male_GutMic
## W = 109, p-value = 0.5676
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                         Sex == "Female", q1, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                          Sex == "Male", q1, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  spi_fem_GutMic and spi_male_GutMic
## W = 55, p-value = 0.8596
## alternative hypothesis: true location shift is not equal to 0

Differences in snout vent length (SVL) among lizard species

set.seed(1234)
Data_SVL <- subset(Micro_div, select = c(Species_Scel, SVL))
kruskal.t <- Data_SVL %>% kruskal_test(SVL ~ Species_Scel)
kruskal.t
Post_hoc  <- Data_SVL %>% dunn_test(SVL ~ Species_Scel, p.adjust.method = "bonferroni") 
Post_hoc

Visualization: box plots with p-values

Post_hoc <- Post_hoc %>% 
  add_xy_position(x = "Species_Scel")
ggboxplot(Data_SVL, x = "Species_Scel", y = "SVL", fill = "Species_Scel") +
  xlab(element_blank())+
  scale_fill_manual(values = c("#56B4E9","#009E73", "#999999","#E69F00"))+
      theme_classic() +
   theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank()) +
  theme(legend.text = element_text(face = "italic"))+
  stat_pvalue_manual(Post_hoc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(kruskal.t, detailed = TRUE),
    caption = get_pwc_label(Post_hoc))

Generalized Linear Model with quasi-Poisson distribution (q = 1)

M1 <- glm(q1 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
           family = quasipoisson(link = "log"),
           data = Micro_div)
summary(M1)
## 
## Call:
## glm(formula = q1 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                4.262828   0.228362  18.667  < 2e-16
## SEASONRainy                                0.405260   0.151684   2.672 0.009061
## SPECIESSceloporus bicanthalis              0.989264   0.459589   2.152 0.034227
## SPECIESSceloporus grammicus                0.028572   0.168483   0.170 0.865747
## SPECIESSceloporus spinosus                 0.044405   0.217178   0.204 0.838486
## ELEVATION                                 -0.387007   0.399665  -0.968 0.335661
## SVL                                        0.002237   0.003986   0.561 0.576205
## SEQDEPTH                                   0.235518   0.065659   3.587 0.000561
## SEASONRainy:SPECIESSceloporus bicanthalis -0.641160   0.190178  -3.371 0.001132
## SEASONRainy:SPECIESSceloporus grammicus   -0.249934   0.196457  -1.272 0.206811
## SEASONRainy:SPECIESSceloporus spinosus    -0.239895   0.202098  -1.187 0.238567
##                                              
## (Intercept)                               ***
## SEASONRainy                               ** 
## SPECIESSceloporus bicanthalis             *  
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## ELEVATION                                    
## SVL                                          
## SEQDEPTH                                  ***
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.08586)
## 
##     Null deviance: 1537.22  on 94  degrees of freedom
## Residual deviance:  999.94  on 84  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(M1)

M2 <- update(M1,. ~ . -SPECIES:SEASON)
anova(M1, M2, test = "F")
summary(M2)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.462617   0.218130  20.458  < 2e-16 ***
## SEASONRainy                    0.114355   0.078245   1.462  0.14748    
## SPECIESSceloporus bicanthalis  0.652287   0.476593   1.369  0.17463    
## SPECIESSceloporus grammicus   -0.144594   0.116591  -1.240  0.21824    
## SPECIESSceloporus spinosus    -0.114180   0.186455  -0.612  0.54189    
## ELEVATION                     -0.467503   0.427006  -1.095  0.27661    
## SVL                            0.001759   0.004152   0.424  0.67284    
## SEQDEPTH                       0.205901   0.068543   3.004  0.00348 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.24543)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1136.0  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
M3 <- update(M2,. ~ . -SVL)
anova(M2, M3, test = "F")
summary(M3)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.53327    0.14023  32.327  < 2e-16 ***
## SEASONRainy                    0.11649    0.07772   1.499  0.13747    
## SPECIESSceloporus bicanthalis  0.66616    0.47539   1.401  0.16464    
## SPECIESSceloporus grammicus   -0.11718    0.09653  -1.214  0.22803    
## SPECIESSceloporus spinosus    -0.04897    0.10416  -0.470  0.63943    
## ELEVATION                     -0.48877    0.42406  -1.153  0.25220    
## SEQDEPTH                       0.20354    0.06796   2.995  0.00356 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.13221)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1138.2  on 88  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
M4 <- update(M3,. ~ . -ELEVATION) 
anova(M3, M4, test = "F")
summary(M4)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.66137    0.08481  54.964   <2e-16 ***
## SEASONRainy                    0.11993    0.07821   1.533   0.1287    
## SPECIESSceloporus bicanthalis  0.12860    0.09341   1.377   0.1721    
## SPECIESSceloporus grammicus   -0.11600    0.09699  -1.196   0.2349    
## SPECIESSceloporus spinosus    -0.01756    0.10109  -0.174   0.8625    
## SEQDEPTH                       0.20435    0.06833   2.991   0.0036 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.2577)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1157.3  on 89  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
confint(M4)
## Waiting for profiling to be done...
##                                     2.5 %     97.5 %
## (Intercept)                    4.49221597 4.82475124
## SEASONRainy                   -0.03336187 0.27330751
## SPECIESSceloporus bicanthalis -0.05390201 0.31251946
## SPECIESSceloporus grammicus   -0.30596309 0.07452011
## SPECIESSceloporus spinosus    -0.21612307 0.18046386
## SEQDEPTH                       0.06781350 0.33576942
M5 <- update(M4,. ~ . -SEQDEPTH) 
anova(M4, M5, test="F")
summary(M5)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.60321    0.08660  53.155  < 2e-16 ***
## SEASONRainy                    0.22149    0.07245   3.057  0.00294 ** 
## SPECIESSceloporus bicanthalis  0.15745    0.09734   1.618  0.10927    
## SPECIESSceloporus grammicus   -0.11983    0.10129  -1.183  0.23993    
## SPECIESSceloporus spinosus    -0.01096    0.10574  -0.104  0.91766    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 13.41224)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1260.8  on 90  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Best Model
M6 <- glm(q1 ~ SEASON*SPECIES,
          family = quasipoisson(link = "log"),
          data = Micro_div)
anova(M6, test="F")
summary(M6)
## 
## Call:
## glm(formula = q1 ~ SEASON * SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 4.3930     0.1316  33.371  < 2e-16
## SEASONRainy                                 0.5318     0.1561   3.406  0.00100
## SPECIESSceloporus bicanthalis               0.5414     0.1656   3.270  0.00154
## SPECIESSceloporus grammicus                 0.1005     0.1635   0.615  0.54032
## SPECIESSceloporus spinosus                  0.1593     0.1715   0.929  0.35550
## SEASONRainy:SPECIESSceloporus bicanthalis  -0.5796     0.2020  -2.869  0.00517
## SEASONRainy:SPECIESSceloporus grammicus    -0.3298     0.2065  -1.597  0.11392
## SEASONRainy:SPECIESSceloporus spinosus     -0.2364     0.2153  -1.098  0.27537
##                                              
## (Intercept)                               ***
## SEASONRainy                               ** 
## SPECIESSceloporus bicanthalis             ** 
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.61408)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1151.1  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
confint(M6)
## Waiting for profiling to be done...
##                                                2.5 %      97.5 %
## (Intercept)                                4.1233396  4.64031820
## SEASONRainy                                0.2317199  0.84503460
## SPECIESSceloporus bicanthalis              0.2207910  0.87128269
## SPECIESSceloporus grammicus               -0.2155594  0.42677045
## SPECIESSceloporus spinosus                -0.1742816  0.49967406
## SEASONRainy:SPECIESSceloporus bicanthalis -0.9789785 -0.18613994
## SEASONRainy:SPECIESSceloporus grammicus   -0.7385988  0.07179439
## SEASONRainy:SPECIESSceloporus spinosus    -0.6616867  0.18335432
coef(M6)
##                               (Intercept) 
##                                 4.3929537 
##                               SEASONRainy 
##                                 0.5317728 
##             SPECIESSceloporus bicanthalis 
##                                 0.5413804 
##               SPECIESSceloporus grammicus 
##                                 0.1005116 
##                SPECIESSceloporus spinosus 
##                                 0.1593282 
## SEASONRainy:SPECIESSceloporus bicanthalis 
##                                -0.5795744 
##   SEASONRainy:SPECIESSceloporus grammicus 
##                                -0.3297727 
##    SEASONRainy:SPECIESSceloporus spinosus 
##                                -0.2363617

Residual plot of the best model

Plot model assumption

plot_model(M6, type = "eff", terms = "SPECIES")

plot_model(M6, show.values = TRUE, value.offset = .3, width = 0.1,
           vline.color = "#E69F00")

residuals_q1 <- plot_model(M6, type = "eff", terms = c("SEASON","SPECIES"),
                           colors = c("#999999", "#E69F00", "#56B4E9", 
                                      "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))
residuals_q1 

ggsave("../figures/residuals_q1.jpeg", width = 5.0, height = 3.0, dpi = 300)

Contrasts (dry vs rainy), CI95%, p values

Tukey Test

library(emmeans)
emmeans(M6, specs = ~ SPECIES*SEASON) %>%
  contrast() %>%
  as.data.frame()
emm <- emmeans(M6, spec = ~ SPECIES*SEASON, 
               type = "response")
print(emm)
##  SPECIES                SEASON  rate    SE  df asymp.LCL asymp.UCL
##  Sceloporus aeneus      Dry     80.9 10.65 Inf      62.5       105
##  Sceloporus bicanthalis Dry    139.0 13.96 Inf     114.1       169
##  Sceloporus grammicus   Dry     89.4  8.67 Inf      74.0       108
##  Sceloporus spinosus    Dry     94.8 10.43 Inf      76.5       118
##  Sceloporus aeneus      Rainy  137.7 11.56 Inf     116.8       162
##  Sceloporus bicanthalis Rainy  132.5 10.56 Inf     113.3       155
##  Sceloporus grammicus   Rainy  109.4 10.31 Inf      91.0       132
##  Sceloporus spinosus    Rainy  127.4 12.68 Inf     104.9       155
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
X <- contrast(emm, method = "pairwise")
confint(X)
M1.emm.s <- emmeans(M6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s, adjust = "tukey", infer=c(TRUE,TRUE))
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry        -0.54138 0.166 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry          -0.10051 0.164 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry           -0.15933 0.172 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy           -0.53177 0.156 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy      -0.49358 0.154 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy        -0.30251 0.162 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy         -0.45474 0.165 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry      0.44087 0.140 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry       0.38205 0.149 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy       0.00961 0.131 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy  0.04780 0.128 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy    0.23887 0.138 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy     0.08664 0.141 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry        -0.05882 0.147 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy        -0.43126 0.128 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy   -0.39307 0.125 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy     -0.20200 0.135 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy      -0.35423 0.139 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy         -0.37244 0.138 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy    -0.33425 0.136 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy      -0.14318 0.145 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy       -0.29541 0.148 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy     0.03819 0.116 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy       0.22926 0.126 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy        0.07703 0.130 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy  0.19107 0.123 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy   0.03884 0.127 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy    -0.15223 0.137 Inf
##  asymp.LCL asymp.UCL z.ratio p.value
##    -1.0432   -0.0396  -3.270  0.0239
##    -0.5961    0.3950  -0.615  0.9987
##    -0.6792    0.3605  -0.929  0.9833
##    -1.0050   -0.0585  -3.406  0.0152
##    -0.9599   -0.0272  -3.208  0.0292
##    -0.7931    0.1880  -1.869  0.5723
##    -0.9548    0.0454  -2.756  0.1064
##     0.0178    0.8640   3.158  0.0341
##    -0.0693    0.8334   2.566  0.1686
##    -0.3871    0.4063   0.073  1.0000
##    -0.3407    0.4363   0.373  1.0000
##    -0.1784    0.6561   1.735  0.6640
##    -0.3418    0.5151   0.613  0.9987
##    -0.5032    0.3855  -0.401  0.9999
##    -0.8200   -0.0425  -3.362  0.0176
##    -0.7734   -0.0127  -3.132  0.0369
##    -0.6117    0.2077  -1.495  0.8107
##    -0.7753    0.0668  -2.550  0.1748
##    -0.7917    0.0469  -2.692  0.1248
##    -0.7458    0.0773  -2.462  0.2121
##    -0.5819    0.2956  -0.989  0.9761
##    -0.7448    0.1540  -1.992  0.4871
##    -0.3126    0.3890   0.330  1.0000
##    -0.1531    0.6116   1.817  0.6082
##    -0.3175    0.4716   0.592  0.9990
##    -0.1828    0.5649   1.549  0.7805
##    -0.3475    0.4251   0.305  1.0000
##    -0.5674    0.2629  -1.111  0.9546
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 8 estimates 
## P value adjustment: tukey method for comparing a family of 8 estimates
pwpm(M1.emm.s)
##                              Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry                       [4.39]                     0.0239
## Sceloporus bicanthalis Dry                -0.54138                     [4.93]
## Sceloporus grammicus Dry                  -0.10051                    0.44087
## Sceloporus spinosus Dry                   -0.15933                    0.38205
## Sceloporus aeneus Rainy                   -0.53177                    0.00961
## Sceloporus bicanthalis Rainy              -0.49358                    0.04780
## Sceloporus grammicus Rainy                -0.30251                    0.23887
## Sceloporus spinosus Rainy                 -0.45474                    0.08664
##                              Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry                          0.9987                  0.9833
## Sceloporus bicanthalis Dry                     0.0341                  0.1686
## Sceloporus grammicus Dry                       [4.49]                  0.9999
## Sceloporus spinosus Dry                      -0.05882                  [4.55]
## Sceloporus aeneus Rainy                      -0.43126                -0.37244
## Sceloporus bicanthalis Rainy                 -0.39307                -0.33425
## Sceloporus grammicus Rainy                   -0.20200                -0.14318
## Sceloporus spinosus Rainy                    -0.35423                -0.29541
##                              Sceloporus aeneus Rainy
## Sceloporus aeneus Dry                         0.0152
## Sceloporus bicanthalis Dry                    1.0000
## Sceloporus grammicus Dry                      0.0176
## Sceloporus spinosus Dry                       0.1248
## Sceloporus aeneus Rainy                       [4.92]
## Sceloporus bicanthalis Rainy                 0.03819
## Sceloporus grammicus Rainy                   0.22926
## Sceloporus spinosus Rainy                    0.07703
##                              Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry                              0.0292
## Sceloporus bicanthalis Dry                         1.0000
## Sceloporus grammicus Dry                           0.0369
## Sceloporus spinosus Dry                            0.2121
## Sceloporus aeneus Rainy                            1.0000
## Sceloporus bicanthalis Rainy                       [4.89]
## Sceloporus grammicus Rainy                        0.19107
## Sceloporus spinosus Rainy                         0.03884
##                              Sceloporus grammicus Rainy
## Sceloporus aeneus Dry                            0.5723
## Sceloporus bicanthalis Dry                       0.6640
## Sceloporus grammicus Dry                         0.8107
## Sceloporus spinosus Dry                          0.9761
## Sceloporus aeneus Rainy                          0.6082
## Sceloporus bicanthalis Rainy                     0.7805
## Sceloporus grammicus Rainy                       [4.70]
## Sceloporus spinosus Rainy                      -0.15223
##                              Sceloporus spinosus Rainy
## Sceloporus aeneus Dry                           0.1064
## Sceloporus bicanthalis Dry                      0.9987
## Sceloporus grammicus Dry                        0.1748
## Sceloporus spinosus Dry                         0.4871
## Sceloporus aeneus Rainy                         0.9990
## Sceloporus bicanthalis Rainy                    1.0000
## Sceloporus grammicus Rainy                      0.9546
## Sceloporus spinosus Rainy                       [4.85]
## 
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values   adjust = "tukey"
## Diagonal: [Estimates] (emmean) 
## Lower triangle: Comparisons (estimate)   earlier vs. later
summary(M1.emm.s, type="response")
plot(M1.emm.s, comparisons=TRUE, type="response")

Evaluate whether sex variable influence on gut microbiota diversity (q=2)

# Paired test (Wilcoxon test, q2) 

# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                            Sex == "Female", q2, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                             Sex == "Male", q2, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  aeneus_fem_GutMic and aeneus_male_GutMic
## W = 39, p-value = 0.1713
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                          Sex == "Female", q2, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                           Sex == "Male", q2, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                          Sex == "Female", q2, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                           Sex == "Male", q2, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  gram_fem_GutMic and gram_male_GutMic
## W = 107, p-value = 0.6313
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                         Sex == "Female", q2, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                          Sex == "Male", q2, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  spi_fem_GutMic and spi_male_GutMic
## W = 40, p-value = 0.4137
## alternative hypothesis: true location shift is not equal to 0

Generalized Linear Model with quasi-Poisson distribution (q = 2)

MODEL1 <- glm(q2 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
            family = quasipoisson(link = "log"),
            data = Micro_div)
summary(MODEL1)
## 
## Call:
## glm(formula = q2 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                3.489e+00  3.263e-01  10.694
## SEASONRainy                                3.432e-01  2.143e-01   1.601
## SPECIESSceloporus bicanthalis              1.267e+00  6.553e-01   1.933
## SPECIESSceloporus grammicus                6.931e-02  2.358e-01   0.294
## SPECIESSceloporus spinosus                 1.270e-01  3.062e-01   0.415
## ELEVATION                                 -4.765e-01  5.739e-01  -0.830
## SVL                                       -2.387e-05  5.754e-03  -0.004
## SEQDEPTH                                   1.204e-01  1.009e-01   1.193
## SEASONRainy:SPECIESSceloporus bicanthalis -8.202e-01  2.652e-01  -3.093
## SEASONRainy:SPECIESSceloporus grammicus   -1.219e-01  2.754e-01  -0.442
## SEASONRainy:SPECIESSceloporus spinosus    -3.317e-01  2.909e-01  -1.140
##                                           Pr(>|t|)    
## (Intercept)                                < 2e-16 ***
## SEASONRainy                                0.11305    
## SPECIESSceloporus bicanthalis              0.05657 .  
## SPECIESSceloporus grammicus                0.76956    
## SPECIESSceloporus spinosus                 0.67944    
## ELEVATION                                  0.40867    
## SVL                                        0.99670    
## SEQDEPTH                                   0.23636    
## SEASONRainy:SPECIESSceloporus bicanthalis  0.00269 ** 
## SEASONRainy:SPECIESSceloporus grammicus    0.65931    
## SEASONRainy:SPECIESSceloporus spinosus     0.25746    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 9.481509)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  818.79  on 84  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(MODEL1)

MODEL2 <- update(MODEL1,. ~ . -SPECIES:SEASON)
anova(MODEL1, MODEL2, test = "F")
summary(MODEL2)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.708501   0.314372  11.797   <2e-16 ***
## SEASONRainy                    0.025511   0.109991   0.232    0.817    
## SPECIESSceloporus bicanthalis  0.933649   0.676327   1.380    0.171    
## SPECIESSceloporus grammicus   -0.025499   0.165697  -0.154    0.878    
## SPECIESSceloporus spinosus    -0.051587   0.267991  -0.192    0.848    
## ELEVATION                     -0.626981   0.606572  -1.034    0.304    
## SVL                           -0.001225   0.005989  -0.205    0.838    
## SEQDEPTH                       0.041928   0.106778   0.393    0.696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.47637)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  934.33  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL3 <- update(MODEL2,. ~ . -SVL)
anova(MODEL2, MODEL3, test="F")
summary(MODEL3)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.65897    0.19913  18.375   <2e-16 ***
## SEASONRainy                    0.02422    0.10921   0.222    0.825    
## SPECIESSceloporus bicanthalis  0.92494    0.66933   1.382    0.171    
## SPECIESSceloporus grammicus   -0.04438    0.13685  -0.324    0.746    
## SPECIESSceloporus spinosus    -0.09673    0.15179  -0.637    0.526    
## ELEVATION                     -0.61281    0.59747  -1.026    0.308    
## SEQDEPTH                       0.04319    0.10604   0.407    0.685    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.36163)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  934.77  on 88  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL4 <- update(MODEL3,. ~ . -SEQDEPTH)
anova(MODEL3, MODEL4, test="F")
summary(MODEL4)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.64561    0.19533  18.664   <2e-16 ***
## SEASONRainy                    0.04488    0.09589   0.468    0.641    
## SPECIESSceloporus bicanthalis  0.93231    0.66505   1.402    0.164    
## SPECIESSceloporus grammicus   -0.04395    0.13615  -0.323    0.748    
## SPECIESSceloporus spinosus    -0.09462    0.15095  -0.627    0.532    
## ELEVATION                     -0.61383    0.59391  -1.034    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.26101)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  936.46  on 89  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL5 <- update(MODEL4,. ~ . -ELEVATION) # Best model
anova(MODEL4, MODEL5, test="F")
summary(MODEL5)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.80662    0.11637  32.712   <2e-16 ***
## SEASONRainy                    0.04940    0.09642   0.512   0.6097    
## SPECIESSceloporus bicanthalis  0.25823    0.13175   1.960   0.0531 .  
## SPECIESSceloporus grammicus   -0.04288    0.13674  -0.314   0.7546    
## SPECIESSceloporus spinosus    -0.05584    0.14700  -0.380   0.7049    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.3626)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  950.07  on 90  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Best Model
MODEL6 <- glm(q2 ~ SEASON*SPECIES,
              family = quasipoisson(link = "log"),
              data = Micro_div)
summary(MODEL6)
## 
## Call:
## glm(formula = q2 ~ SEASON * SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                3.57869    0.17122  20.902  < 2e-16
## SEASONRainy                                0.40334    0.20706   1.948 0.054645
## SPECIESSceloporus bicanthalis              0.74146    0.20804   3.564 0.000596
## SPECIESSceloporus grammicus                0.08107    0.21339   0.380 0.704917
## SPECIESSceloporus spinosus                 0.15897    0.22310   0.713 0.478026
## SEASONRainy:SPECIESSceloporus bicanthalis -0.79089    0.26304  -3.007 0.003451
## SEASONRainy:SPECIESSceloporus grammicus   -0.14701    0.27125  -0.542 0.589223
## SEASONRainy:SPECIESSceloporus spinosus    -0.32729    0.29016  -1.128 0.262430
##                                              
## (Intercept)                               ***
## SEASONRainy                               .  
## SPECIESSceloporus bicanthalis             ***
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 9.452614)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  840.53  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
coef(MODEL6)
##                               (Intercept) 
##                                3.57869413 
##                               SEASONRainy 
##                                0.40334451 
##             SPECIESSceloporus bicanthalis 
##                                0.74145835 
##               SPECIESSceloporus grammicus 
##                                0.08107416 
##                SPECIESSceloporus spinosus 
##                                0.15897437 
## SEASONRainy:SPECIESSceloporus bicanthalis 
##                               -0.79088870 
##   SEASONRainy:SPECIESSceloporus grammicus 
##                               -0.14700948 
##    SEASONRainy:SPECIESSceloporus spinosus 
##                               -0.32729337
confint(MODEL6)
## Waiting for profiling to be done...
##                                                  2.5 %     97.5 %
## (Intercept)                                3.223230360  3.8964819
## SEASONRainy                                0.006183591  0.8207875
## SPECIESSceloporus bicanthalis              0.342095601  1.1605800
## SPECIESSceloporus grammicus               -0.330307974  0.5093207
## SPECIESSceloporus spinosus                -0.274369614  0.6037278
## SEASONRainy:SPECIESSceloporus bicanthalis -1.313319633 -0.2802605
## SEASONRainy:SPECIESSceloporus grammicus   -0.685102597  0.3801804
## SEASONRainy:SPECIESSceloporus spinosus    -0.902353379  0.2372010

Residual plot of the best model

Plot model assumption

plot_model(MODEL6, type = "eff", terms = "SPECIES")

plot_model(MODEL6, show.values = TRUE, value.offset = .3, width = 0.1,
           vline.color = "#E69F00")

residuals_q2 <- plot_model(MODEL6, type = "eff", terms = c("SEASON","SPECIES"),
                           colors = c("#999999", "#E69F00", "#56B4E9", 
                                      "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))
residuals_q2 

ggsave("../figures/residuals_q2.jpeg", width = 5.0, height = 3.0, dpi = 300)

Contrasts (dry vs rainy), CI95%, p values

## Tukey Test

emmeans(MODEL6, specs = ~ SPECIES*SEASON) %>%
  contrast() %>%
  as.data.frame()
emm_q2 <- emmeans(MODEL6, spec = ~ SPECIES*SEASON, 
                  type = "response")
print(emm_q2)
##  SPECIES                SEASON rate   SE  df asymp.LCL asymp.UCL
##  Sceloporus aeneus      Dry    35.8 6.13 Inf      25.6      50.1
##  Sceloporus bicanthalis Dry    75.2 8.89 Inf      59.7      94.8
##  Sceloporus grammicus   Dry    38.9 4.95 Inf      30.3      49.9
##  Sceloporus spinosus    Dry    42.0 6.01 Inf      31.7      55.6
##  Sceloporus aeneus      Rainy  53.6 6.24 Inf      42.7      67.4
##  Sceloporus bicanthalis Rainy  51.0 5.67 Inf      41.1      63.5
##  Sceloporus grammicus   Rainy  50.2 6.04 Inf      39.7      63.6
##  Sceloporus spinosus    Rainy  45.3 6.55 Inf      34.1      60.1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
X_q2 <- contrast(emm_q2, method = "pairwise")
confint(X_q2)
M1.emm.s_q2 <- emmeans(MODEL6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s_q2, adjust = "tukey", infer=c(TRUE,TRUE))
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry         -0.7415 0.208 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry           -0.0811 0.213 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry            -0.1590 0.223 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy            -0.4033 0.207 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy       -0.3539 0.204 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy         -0.3374 0.209 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy          -0.2350 0.224 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry       0.6604 0.174 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry        0.5825 0.186 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy        0.3381 0.166 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy   0.3875 0.162 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy     0.4040 0.169 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy      0.5064 0.187 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry         -0.0779 0.192 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy         -0.3223 0.173 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy    -0.2728 0.169 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy      -0.2563 0.175 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy       -0.1540 0.193 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy          -0.2444 0.184 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy     -0.1949 0.181 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy       -0.1784 0.187 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy        -0.0761 0.203 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy      0.0494 0.161 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy        0.0659 0.167 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy         0.1683 0.186 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy   0.0165 0.164 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy    0.1189 0.182 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy      0.1024 0.188 Inf
##  asymp.LCL asymp.UCL z.ratio p.value
##    -1.3720    -0.111  -3.564  0.0087
##    -0.7278     0.566  -0.380  0.9999
##    -0.8352     0.517  -0.713  0.9967
##    -1.0309     0.224  -1.948  0.5176
##    -0.9726     0.265  -1.734  0.6649
##    -0.9717     0.297  -1.612  0.7431
##    -0.9139     0.444  -1.049  0.9668
##     0.1338     1.187   3.801  0.0036
##     0.0201     1.145   3.139  0.0361
##    -0.1647     0.841   2.038  0.4560
##    -0.1041     0.879   2.389  0.2465
##    -0.1072     0.915   2.396  0.2433
##    -0.0592     1.072   2.714  0.1183
##    -0.6584     0.503  -0.407  0.9999
##    -0.8453     0.201  -1.868  0.5734
##    -0.7851     0.239  -1.614  0.7418
##    -0.7874     0.275  -1.463  0.8272
##    -0.7376     0.430  -0.800  0.9932
##    -0.8034     0.315  -1.325  0.8897
##    -0.7439     0.354  -1.076  0.9618
##    -0.7450     0.388  -0.955  0.9805
##    -0.6921     0.540  -0.374  1.0000
##    -0.4384     0.537   0.307  1.0000
##    -0.4416     0.573   0.394  0.9999
##    -0.3940     0.731   0.907  0.9855
##    -0.4799     0.513   0.101  1.0000
##    -0.4334     0.671   0.652  0.9981
##    -0.4674     0.672   0.545  0.9994
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 8 estimates 
## P value adjustment: tukey method for comparing a family of 8 estimates
pairs(M1.emm.s_q2, adjust = "tukey")
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry         -0.7415 0.208 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry           -0.0811 0.213 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry            -0.1590 0.223 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy            -0.4033 0.207 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy       -0.3539 0.204 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy         -0.3374 0.209 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy          -0.2350 0.224 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry       0.6604 0.174 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry        0.5825 0.186 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy        0.3381 0.166 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy   0.3875 0.162 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy     0.4040 0.169 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy      0.5064 0.187 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry         -0.0779 0.192 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy         -0.3223 0.173 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy    -0.2728 0.169 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy      -0.2563 0.175 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy       -0.1540 0.193 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy          -0.2444 0.184 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy     -0.1949 0.181 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy       -0.1784 0.187 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy        -0.0761 0.203 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy      0.0494 0.161 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy        0.0659 0.167 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy         0.1683 0.186 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy   0.0165 0.164 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy    0.1189 0.182 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy      0.1024 0.188 Inf
##  z.ratio p.value
##   -3.564  0.0087
##   -0.380  0.9999
##   -0.713  0.9967
##   -1.948  0.5176
##   -1.734  0.6649
##   -1.612  0.7431
##   -1.049  0.9668
##    3.801  0.0036
##    3.139  0.0361
##    2.038  0.4560
##    2.389  0.2465
##    2.396  0.2433
##    2.714  0.1183
##   -0.407  0.9999
##   -1.868  0.5734
##   -1.614  0.7418
##   -1.463  0.8272
##   -0.800  0.9932
##   -1.325  0.8897
##   -1.076  0.9618
##   -0.955  0.9805
##   -0.374  1.0000
##    0.307  1.0000
##    0.394  0.9999
##    0.907  0.9855
##    0.101  1.0000
##    0.652  0.9981
##    0.545  0.9994
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 8 estimates
pwpm(M1.emm.s_q2)
##                              Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry                       [3.58]                     0.0087
## Sceloporus bicanthalis Dry                 -0.7415                     [4.32]
## Sceloporus grammicus Dry                   -0.0811                     0.6604
## Sceloporus spinosus Dry                    -0.1590                     0.5825
## Sceloporus aeneus Rainy                    -0.4033                     0.3381
## Sceloporus bicanthalis Rainy               -0.3539                     0.3875
## Sceloporus grammicus Rainy                 -0.3374                     0.4040
## Sceloporus spinosus Rainy                  -0.2350                     0.5064
##                              Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry                          0.9999                  0.9967
## Sceloporus bicanthalis Dry                     0.0036                  0.0361
## Sceloporus grammicus Dry                       [3.66]                  0.9999
## Sceloporus spinosus Dry                       -0.0779                  [3.74]
## Sceloporus aeneus Rainy                       -0.3223                 -0.2444
## Sceloporus bicanthalis Rainy                  -0.2728                 -0.1949
## Sceloporus grammicus Rainy                    -0.2563                 -0.1784
## Sceloporus spinosus Rainy                     -0.1540                 -0.0761
##                              Sceloporus aeneus Rainy
## Sceloporus aeneus Dry                         0.5176
## Sceloporus bicanthalis Dry                    0.4560
## Sceloporus grammicus Dry                      0.5734
## Sceloporus spinosus Dry                       0.8897
## Sceloporus aeneus Rainy                       [3.98]
## Sceloporus bicanthalis Rainy                  0.0494
## Sceloporus grammicus Rainy                    0.0659
## Sceloporus spinosus Rainy                     0.1683
##                              Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry                              0.6649
## Sceloporus bicanthalis Dry                         0.2465
## Sceloporus grammicus Dry                           0.7418
## Sceloporus spinosus Dry                            0.9618
## Sceloporus aeneus Rainy                            1.0000
## Sceloporus bicanthalis Rainy                       [3.93]
## Sceloporus grammicus Rainy                         0.0165
## Sceloporus spinosus Rainy                          0.1189
##                              Sceloporus grammicus Rainy
## Sceloporus aeneus Dry                            0.7431
## Sceloporus bicanthalis Dry                       0.2433
## Sceloporus grammicus Dry                         0.8272
## Sceloporus spinosus Dry                          0.9805
## Sceloporus aeneus Rainy                          0.9999
## Sceloporus bicanthalis Rainy                     1.0000
## Sceloporus grammicus Rainy                       [3.92]
## Sceloporus spinosus Rainy                        0.1024
##                              Sceloporus spinosus Rainy
## Sceloporus aeneus Dry                           0.9668
## Sceloporus bicanthalis Dry                      0.1183
## Sceloporus grammicus Dry                        0.9932
## Sceloporus spinosus Dry                         1.0000
## Sceloporus aeneus Rainy                         0.9855
## Sceloporus bicanthalis Rainy                    0.9981
## Sceloporus grammicus Rainy                      0.9994
## Sceloporus spinosus Rainy                       [3.81]
## 
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values   adjust = "tukey"
## Diagonal: [Estimates] (emmean) 
## Lower triangle: Comparisons (estimate)   earlier vs. later
summary(M1.emm.s_q2, type="response")
plot(M1.emm.s_q2, comparisons=TRUE, type="response")

Relation between alpha taxonomic microbiota and diet

Micro_div$Ind <- substr(Micro_div$SampleID, 8, nchar(Micro_div$SampleID) - 0)

diet_phylo_div <- read.csv("../data/diet_Hill_numbers_q012.csv", header = TRUE) 
# Generar una nueva columna con la información de otra columna, excluyendo los últimos 4 caracteres
diet_phylo_div$Ind <- substr(diet_phylo_div$SampleID, 8, nchar(diet_phylo_div$SampleID) - 4)
#names(diet_phylo_div)[names(diet_phylo_div) %in% c("q0", "q1", "q2")] <- c("d.p.0", "d.p.1", #"d.p.2")
divs <- Micro_div %>% 
  inner_join(diet_phylo_div, by = c("Ind"="Ind"))

#write.table(divs, file="../data/full_divs.txt", sep = "\t")

##Plot

taxq1_tq0 <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), xlab= "Diet taxonomic diversity", ylab="Gut bacterial diversity (q=1)")+
  theme(legend.title = element_blank(), legend.position = "bottom")


taxq2_tq0 <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Diet taxonomic diversity", ylab="Gut bacterial diversity (q=2)")+
  theme(legend.title = element_blank(), legend.position = "bottom")


phyq1_pq0 <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") , 
   xlab= " Diet phylogenetic diversity", ylab="Gut bacterial diversity (q=1)")+
  theme(legend.title = element_blank(), legend.position = "bottom")


phyq2_pq0 <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Diet phylogenetic diversity", ylab="Gut bacterial diversity (q=2)")+
  theme(legend.title = element_blank(), legend.position = "bottom")


taxq1_tq0.sp <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Diet taxonomic diversity", ylab="Gut bacterial diversity (q=1)",
   facet.by = "Species_Scel" #,
   #add = "reg.line",  # Add regressin line
   #add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   #conf.int = TRUE, # Add confidence interval
   #cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   #cor.coeff.args = list(method = "spearman", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none")


taxq2_tq0.sp <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Diet taxonomic diversity", ylab="Gut bacterial diversity (q=2)",
   facet.by = "Species_Scel" #,
   #add = "reg.line",  # Add regressin line
   #add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   #conf.int = TRUE, # Add confidence interval
   #cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   #cor.coeff.args = list(method = "spearman", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none")


phyq1_pq0.sp <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") , 
   xlab= " Diet phylogenetic diversity", ylab="Gut bacterial diversity (q=1)",
   facet.by = "Species_Scel"#,
   #add = "reg.line",  # Add regressin line
   #add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   #conf.int = TRUE, # Add confidence interval
   #cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   #cor.coeff.args = list(method = "spearman", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none")


phyq2_pq0.sp <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Diet phylogenetic diversity", ylab="Gut bacterial diversity (q=2)", 
   facet.by = "Species_Scel" #,
   #add = "reg.line",  # Add regressin line
   #add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   #conf.int = TRUE, # Add confidence interval
   #cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   #cor.coeff.args = list(method = "spearman", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none")

taxq1_tq0

taxq1_tq0.sp

taxq2_tq0

taxq2_tq0.sp

phyq1_pq0

phyq1_pq0.sp

phyq2_pq0

phyq2_pq0.sp

diet.vs.micro1<- plot_grid(taxq1_tq0, taxq1_tq0.sp, 
                          labels = "AUTO", ncol = 1)
ggsave("../figures/dietvsmicro1.jpeg", width=8, height=9, dpi=300)
diet.vs.micro2<- plot_grid(taxq2_tq0, taxq2_tq0.sp, 
                          labels = "AUTO", ncol = 1)
ggsave("../figures/dietvsmicro2.jpeg", width=8, height=9, dpi=300)
diet.vs.micro3<- plot_grid(phyq1_pq0, phyq1_pq0.sp, 
                          labels = "AUTO", ncol = 1)
ggsave("../figures/dietvsmicro3.jpeg", width=8, height=9, dpi=300)
diet.vs.micro4<- plot_grid(phyq2_pq0, phyq2_pq0.sp, 
                          labels = "AUTO", ncol = 1)
ggsave("../figures/dietvsmicro4.jpeg", width=8, height=9, dpi=300)

Effect of diet on gut microbiota

loading files

richness_microb <- read.csv("../data/Hill_numbers_q012.csv", header = TRUE) %>% 
  dplyr::select(SampleID, q0, q1, q2)
metadata <- read.csv("../data/Metadata_model_DietMic.csv", header = TRUE, check.names = F)
Microb_data <- richness_microb %>% inner_join(metadata, by = c("SampleID"="SampleID"))

Check distribution of the dependent variable (diversity q=0, Hill numbers)

shapiro.test(Microb_data$q0) # Data are not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  Microb_data$q0
## W = 0.92649, p-value = 0.0009384
shapiro.test(Microb_data$q1) # Normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  Microb_data$q1
## W = 0.97555, p-value = 0.2334
hist(Microb_data$q0, col = "darkgreen") # q0

hist(log(Microb_data$q0), col = "orange") # log transformation

hist(Microb_data$q1, col = "gray") # q1

hist(Microb_data$q2, col = "blue") # q2

bartlett.test(q0 ~ Species_Scel, data = Microb_data) 
## 
##  Bartlett test of homogeneity of variances
## 
## data:  q0 by Species_Scel
## Bartlett's K-squared = 9.0931, df = 3, p-value = 0.02808
leveneTest(q0 ~ Species_Scel, data = Microb_data)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

Plot variances among lizard species (q0)

boxplot(q0 ~ Species_Scel, data = Microb_data,
        main = "Differences in variance among lizard species",
        xlab = "Species",
        ylab = " Hill numbers (q0)",
        col = "#4682B433",
        border = "black")

###Preparing the data

#Declare as factors
SPECIES <- as.factor(Microb_data$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(Microb_data$Season) # Dry and Rainy
SEX <- as.factor(Microb_data$Sex) # Male and Female

#Standardize continous independent variables
SVL <- rescale(Microb_data$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(Microb_data$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(Microb_data$SeqDepth, binary.inputs = "center")

Linear model with Poisson distribution (taxonomic q = 0)

With Microbial Diversity (q1)

## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ Diet*SPECIES+SEASON,
          family = quasipoisson(link = "log"),
          data = Microb_data)
summary(m1)
## 
## Call:
## glm(formula = q1 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         4.6562151  0.1865713  24.957  < 2e-16 ***
## Diet                               -0.0004683  0.0196703  -0.024  0.98109    
## SPECIESSceloporus bicanthalis       0.0595147  0.2216353   0.269  0.78930    
## SPECIESSceloporus grammicus        -0.1948473  0.2453897  -0.794  0.43059    
## SPECIESSceloporus spinosus          0.1751032  0.2349154   0.745  0.45921    
## SEASONRainy                         0.2689010  0.0857259   3.137  0.00274 ** 
## Diet:SPECIESSceloporus bicanthalis  0.0018106  0.0258806   0.070  0.94448    
## Diet:SPECIESSceloporus grammicus   -0.0084125  0.0312962  -0.269  0.78909    
## Diet:SPECIESSceloporus spinosus    -0.0373851  0.0318303  -1.175  0.24525    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.57198)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 593.29  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(m1)

m2 <- update(m1,. ~ . -Diet:SPECIES)
anova(m1, m2, test="F")
summary(m2)
## 
## Call:
## glm(formula = q1 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.704855   0.130358  36.092  < 2e-16 ***
## Diet                          -0.007877   0.010173  -0.774  0.44187    
## SPECIESSceloporus bicanthalis  0.061506   0.106646   0.577  0.56636    
## SPECIESSceloporus grammicus   -0.256866   0.111237  -2.309  0.02452 *  
## SPECIESSceloporus spinosus    -0.047488   0.116892  -0.406  0.68605    
## SEASONRainy                    0.283586   0.080238   3.534  0.00081 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.42645)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 614.21  on 58  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m3 <- update(m2,. ~ . -Diet)
anova(m2, m3, test="F")
summary(m3) # Best Model
## 
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.63657    0.09649  48.053  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.07079    0.10567   0.670 0.505539    
## SPECIESSceloporus grammicus   -0.24165    0.10918  -2.213 0.030761 *  
## SPECIESSceloporus spinosus    -0.02901    0.11425  -0.254 0.800452    
## SEASONRainy                    0.29232    0.07918   3.692 0.000488 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.37193)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 620.50  on 59  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m4 <- update(m3,. ~ . -SPECIES)
anova(m3, m4, test="F")
summary(m4)
## 
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.55181    0.06272  72.574  < 2e-16 ***
## SEASONRainy  0.33820    0.08005   4.225 7.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.18809)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 728.20  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m5 <- update(m3,. ~ . -SEASON)
anova(m3, m5, test="F")
summary(m5)
## 
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.82632    0.08770  55.033   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.09084    0.11575   0.785   0.4356    
## SPECIESSceloporus grammicus   -0.29554    0.11870  -2.490   0.0156 *  
## SPECIESSceloporus spinosus    -0.08290    0.12430  -0.667   0.5074    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.47322)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 764.28  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))

With Microbial Diversity (q2)

model1 <- glm(q2 ~ Diet*SPECIES+SEASON,
              family = quasipoisson(link = "log"),
              data = Microb_data)
summary(model1)
## 
## Call:
## glm(formula = q2 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         3.62441    0.25085  14.449   <2e-16 ***
## Diet                                0.01764    0.02568   0.687   0.4950    
## SPECIESSceloporus bicanthalis       0.34411    0.29428   1.169   0.2473    
## SPECIESSceloporus grammicus         0.08006    0.32076   0.250   0.8038    
## SPECIESSceloporus spinosus          0.45688    0.31449   1.453   0.1520    
## SEASONRainy                         0.11588    0.10904   1.063   0.2925    
## Diet:SPECIESSceloporus bicanthalis -0.01057    0.03311  -0.319   0.7507    
## Diet:SPECIESSceloporus grammicus   -0.03106    0.04019  -0.773   0.4429    
## Diet:SPECIESSceloporus spinosus    -0.07581    0.04204  -1.803   0.0768 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.255196)
## 
##     Null deviance: 528.4  on 63  degrees of freedom
## Residual deviance: 401.7  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(model1)

model2 <- update(model1,. ~ . -Diet:SPECIES)
anova(model1, model2, test="F")
summary(model2)
## 
## Call:
## glm(formula = q2 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.79639    0.17233  22.030   <2e-16 ***
## Diet                          -0.00507    0.01324  -0.383   0.7030    
## SPECIESSceloporus bicanthalis  0.24216    0.14164   1.710   0.0927 .  
## SPECIESSceloporus grammicus   -0.15074    0.14887  -1.013   0.3155    
## SPECIESSceloporus spinosus    -0.02013    0.15835  -0.127   0.8993    
## SEASONRainy                    0.13451    0.10404   1.293   0.2012    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.364622)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 430.43  on 58  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model3 <- update(model2,. ~ . -Diet)
anova(model2, model3, test="F")
summary(model3)
## 
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.752299   0.127845  29.350   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.247615   0.139789   1.771   0.0817 .  
## SPECIESSceloporus grammicus   -0.140309   0.145194  -0.966   0.3378    
## SPECIESSceloporus spinosus    -0.008757   0.154400  -0.057   0.9550    
## SEASONRainy                    0.140365   0.102034   1.376   0.1741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.248433)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 431.52  on 59  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model4 <- update(model3,. ~ . -SEASON)
anova(model3, model4, test="F") # Best model
summary(model4)
## 
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.84098    0.10960  35.046   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.25751    0.13984   1.841   0.0705 .  
## SPECIESSceloporus grammicus   -0.16641    0.14420  -1.154   0.2531    
## SPECIESSceloporus spinosus    -0.03486    0.15350  -0.227   0.8211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.272215)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 445.32  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model5 <- update(model4,. ~ . -SPECIES)
anova(model4, model5, test="F")
summary(model5)
## 
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = Microb_data)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.8565     0.0515   74.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 8.029905)
## 
##     Null deviance: 528.4  on 63  degrees of freedom
## Residual deviance: 528.4  on 63  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

Linear models with Poisson distribution (Phylogenetic q = 0)

With Microbial Diversity (q1)

divs$phyDiet<- divs$pq0
#Declare as factors
SPECIES <- as.factor(divs$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(divs$Season) # Dry and Rainy
SEX <- as.factor(divs$Sex) # Male and Female

#Standardize continous independent variables
SVL <- rescale(divs$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(divs$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(divs$SeqDepth, binary.inputs = "center")

## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ phyDiet*SPECIES+SEASON,
          family = quasipoisson(link = "log"),
          data = divs)
summary(m1)
## 
## Call:
## glm(formula = q1 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            4.661408   0.206567  22.566  < 2e-16 ***
## phyDiet                               -0.005446   0.069609  -0.078  0.93791    
## SPECIESSceloporus bicanthalis         -0.042498   0.236501  -0.180  0.85803    
## SPECIESSceloporus grammicus           -0.262152   0.243544  -1.076  0.28628    
## SPECIESSceloporus spinosus            -0.013347   0.247408  -0.054  0.95717    
## SEASONRainy                            0.276381   0.088922   3.108  0.00294 ** 
## phyDiet:SPECIESSceloporus bicanthalis  0.050503   0.084691   0.596  0.55332    
## phyDiet:SPECIESSceloporus grammicus   -0.004840   0.086901  -0.056  0.95577    
## phyDiet:SPECIESSceloporus spinosus    -0.014879   0.082267  -0.181  0.85711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.85698)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 619.92  on 57  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(m1)

m2 <- update(m1,. ~ . -phyDiet:SPECIES)
anova(m1, m2, test="F")
summary(m2)
## 
## Call:
## glm(formula = q1 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.651537   0.133638  34.807  < 2e-16 ***
## phyDiet                        0.001562   0.026733   0.058  0.95360    
## SPECIESSceloporus bicanthalis  0.073228   0.106615   0.687  0.49483    
## SPECIESSceloporus grammicus   -0.274588   0.113016  -2.430  0.01812 *  
## SPECIESSceloporus spinosus    -0.062852   0.109585  -0.574  0.56842    
## SEASONRainy                    0.264348   0.086185   3.067  0.00324 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.49499)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 633.01  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m3 <- update(m2,. ~ . -phyDiet)
anova(m2, m3, test="F")
summary(m3) # Best Model
## 
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.65694    0.09566  48.684  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.07275    0.10543   0.690  0.49280    
## SPECIESSceloporus grammicus   -0.27534    0.11136  -2.472  0.01622 *  
## SPECIESSceloporus spinosus    -0.06223    0.10818  -0.575  0.56720    
## SEASONRainy                    0.26231    0.07814   3.357  0.00136 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.32373)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 633.05  on 61  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m4 <- update(m3,. ~ . -SPECIES)
anova(m3, m4, test="F")
summary(m4)
## 
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.55181    0.06312  72.115  < 2e-16 ***
## SEASONRainy  0.32085    0.07995   4.013  0.00016 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.33081)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 758.28  on 64  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m5 <- update(m3,. ~ . -SEASON)
anova(m3, m5, test="F")
summary(m5)
## 
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.82632    0.08573  56.296  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.09084    0.11315   0.803  0.42514    
## SPECIESSceloporus grammicus   -0.33142    0.11834  -2.801  0.00679 ** 
## SPECIESSceloporus spinosus    -0.08424    0.11604  -0.726  0.47057    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.91967)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 751.44  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))

With Microbial Diversity (q2)

model1 <- glm(q2 ~ phyDiet*SPECIES+SEASON,
              family = quasipoisson(link = "log"),
              data = divs)
summary(model1)
## 
## Call:
## glm(formula = q2 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            3.878108   0.280118  13.845   <2e-16 ***
## phyDiet                               -0.042399   0.096811  -0.438    0.663    
## SPECIESSceloporus bicanthalis          0.005418   0.319419   0.017    0.987    
## SPECIESSceloporus grammicus           -0.258721   0.327267  -0.791    0.432    
## SPECIESSceloporus spinosus             0.045784   0.335106   0.137    0.892    
## SEASONRainy                            0.112807   0.115877   0.974    0.334    
## phyDiet:SPECIESSceloporus bicanthalis  0.101063   0.114513   0.883    0.381    
## phyDiet:SPECIESSceloporus grammicus    0.034825   0.117544   0.296    0.768    
## phyDiet:SPECIESSceloporus spinosus    -0.026253   0.113567  -0.231    0.818    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.702367)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 447.64  on 57  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(model1)

model2 <- update(model1,. ~ . -phyDiet:SPECIES)
anova(model1, model2, test="F")
summary(model2)
## 
## Call:
## glm(formula = q2 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.82070    0.17848  21.407   <2e-16 ***
## phyDiet                       -0.01355    0.03571  -0.379   0.7058    
## SPECIESSceloporus bicanthalis  0.24687    0.14407   1.714   0.0918 .  
## SPECIESSceloporus grammicus   -0.17655    0.15307  -1.153   0.2533    
## SPECIESSceloporus spinosus    -0.04872    0.15269  -0.319   0.7508    
## SEASONRainy                    0.08802    0.11442   0.769   0.4448    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.670608)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 467.11  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model3 <- update(model2,. ~ . -phyDiet)
anova(model2, model3, test="F")
summary(model3)
## 
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.77434    0.12992  29.052   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.24998    0.14286   1.750   0.0852 .  
## SPECIESSceloporus grammicus   -0.17020    0.15117  -1.126   0.2646    
## SPECIESSceloporus spinosus    -0.05479    0.15090  -0.363   0.7178    
## SEASONRainy                    0.10615    0.10332   1.027   0.3083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.57)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 468.22  on 61  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model4 <- update(model3,. ~ . -SEASON)
anova(model3, model4, test="F") # Best model
summary(model4)
## 
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.84098    0.11157  34.427   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.25751    0.14236   1.809   0.0753 .  
## SPECIESSceloporus grammicus   -0.19303    0.14920  -1.294   0.2006    
## SPECIESSceloporus spinosus    -0.06383    0.15031  -0.425   0.6725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.536168)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 476.25  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model5 <- update(model4,. ~ . -SPECIES)
anova(model4, model5, test="F")
summary(model5)
## 
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = divs)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.84306    0.05214    73.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 8.3741)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 568.65  on 65  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))